function [rgbimageFinal] = findBriefcase(rgbImage, xyzImage, closestPoint,videoFrame, displayImages)
%This function finds the briefcase in the image. It then determines the 
%4 corners of the briefcase and estimates the plane for the briefcase.
global meanVals;
%Process the range values to find the image
imageZ = xyzImage(:,:,3);
%figure(1)
colImageZ = reshape(imageZ,1, 307200);
%hist(colImageZ, 50)
%Display the original rgb image
if displayImages>0
%    subplot(3,1,1); 
   imshow(rgbImage); 
end
%Threshold to find a sub-section of the briefcase
%****************************************************************
%Threshold the range image to identify a sub-section of the briefcase
range = 0.1;
[rangeImage, rowsZ, colsZ] = thresholdRangeImage(imageZ,closestPoint, range);

if displayImages>0
%    subplot(3,1,2);
   imshow(rangeImage); 
end
%******************************************************************

%Find the equation of the plane on the briefcase using the subset of points
%******************************************************************
%Generate points to be used for estimating the plane
numPoints = size(rowsZ,1);
permutedPoints = randperm(numPoints);

%Determine the number of sample points
samplePoints = ceil(numPoints/5);
randomPoints = permutedPoints(1:samplePoints);

%Take a sample of the points
colSample = colsZ(randomPoints);
rowSample = rowsZ(randomPoints);

%store the xyz points for the briefcase
points = zeros(samplePoints,3);
briefCasePlane = zeros(480,640);
for ii=1:samplePoints
  points(ii,:)= xyzImage(rowSample(ii),colSample(ii),:);
  briefCasePlane(rowSample(ii),colSample(ii)) =1; 
end

%get the plane coefficients that best fits the points
%[fitlist,plane] = select_patch(points)
plane = fitplane(points);

%Display the briefcase
if displayImages>0
%     subplot(3,1,3);
    imshow(briefCasePlane);
    xlabel('The points used to determine the plane of the briefcase');
end

%************************************************************

%Grow the sub-section of the briefcase to find the whole briefcase
%*****************************************************************
%Find points that include the entire briefcase
[rowsGrow, colsGrow] = find(imageZ >=(closestPoint-0.3) & imageZ < 0);

%store the xyz points for region growing
for ii=1:size(rowsGrow,1)
  allPoints(ii,:)= xyzImage(rowsGrow(ii),colsGrow(ii),:);
end

briefcaseImageFull = zeros(480,640);
for mm=1:size(allPoints,1)
briefcaseImageFull(rowsGrow(mm), colsGrow(mm)) = 1;
end

if displayImages>0
   imshow(briefcaseImageFull); 
   xlabel('The image containing the whole briefcase');
end

%Return the newlist which is all the points that can be
%added to the plane
[newlist,indices] = getallpoints(plane,points,allPoints,allPoints);

%Display the image containing these points
briefcaseImageFull = zeros(480,640);
for ll=1:size(indices,1)
briefcaseImageFull(rowsGrow(indices(ll)), colsGrow(indices(ll))) = 1;
end

if displayImages>0
%     subplot(3,2,1);
   imshow(briefcaseImageFull); 
   xlabel('The image containing the whole briefcase');
end

%Find the rest of the briefcase
%******************************************************************
%perform element-wise multiplication to extract the briefcase
briefcaseImageR = double(rgbImage(:,:,1)).*briefcaseImageFull;

if displayImages>0
%     subplot(3,2,2);
   imshow(briefcaseImageR); 
   xlabel('The thresholded image containing the whole briefcase');
end

%Remove the hand using histogram thresholding
thehist = dohist(briefcaseImageR(:,:,1),11);

%find the mean intensity of the image
thehist(1,:) = [];
numPixels = sum(thehist(:,1));
weights = thehist./numPixels;
vals =(1:255)';
meanIntensity  = sum(vals.*weights);
meanVals = [meanVals; meanIntensity];
%Threshold the image along the red channel to remove the skin
if meanIntensity < 40
threshImage = briefcaseImageR(:,:,1)<120 & briefcaseImageR(:,:,1)>5;
else
threshImage = briefcaseImageR(:,:,1)<80 & briefcaseImageR(:,:,1)>5;
end
%Perform erosion and dilation
SE = strel('square', 4);
openedImage = imopen(threshImage, SE);


%========Added by george
stats = regionprops(openedImage, 'area');
%disp(size(stats));

if (length(stats)>0)
    for i=1:length(stats)
        areas(i)=stats(i).Area;
    end
    %areas = stats.Area
    [maxArea, I] = max(areas); 
    openedImage = bwareaopen(openedImage,maxArea-2);

end

if displayImages>0
% subplot(3,2,3);
imshow(openedImage);
xlabel('The briefcase after hand removal');
hold on
end

[B, L] = bwboundaries(openedImage,'noholes');
boundary = B{1};

eq = zeros(4,2);
%fign = uint8(zeros(480,640));
for i = 1:4
    
    [flag,t,d,nr,nc,count,frl,fcl,newcountl] = ransacline(boundary(:,1),boundary(:,2),2,0.25,0.01,0.001,60,3);
    %Remove boundary points used to construct the current ransac line
    boundary=setdiff(boundary,[frl fcl],'rows');
    
    %Plot the ransac line
    [cr,cc] = plotline(t,d);
    plot(cc,cr)
    
    %Calculate the equation of the line
    eq(i,:) = [sin(t)/cos(t) d/cos(t)];
    
end


hold off
%Now we have the corner points. Figure out which
%ones correspond to which corner point
stats = regionprops(openedImage, 'Centroid','BoundingBox');
centroid = stats.Centroid;
%centroid = permute(centroid,[2,1]);
centroid = centroid([2 1]);
corners = zeros(4,2);
for i=1:4
    for j=1:4
        
       %compute the intersecting point
       A = [eq(i,1) -1; eq(j,1) -1];
       B = [-eq(i,2) -eq(j,2)]';
       point = abs(linsolve(A,B)); 
       
       %if the distance between the intersection point and the 
       %distance = point(1)-centroid(1)
       if norm(point - centroid') < 100
           if (point(1)-centroid(1))*(point(2)-centroid(2)) > 0
              if point(1) < centroid(1)
                  corners(1,:)=point;
              else
                  corners(3,:)=point;
              end
           else
                if point(1) < centroid(1)
                  corners(2,:)=point;
                else
                  corners(4,:)=point;
                end
           end
       end
       
    end 
end

if displayImages>0
% subplot(3,3,1);
imshow(openedImage);
hold on
plot(corners(1,2),corners(1,1),'b+');
plot(corners(2,2),corners(2,1),'b+');
plot(corners(3,2),corners(3,1),'b+');
plot(corners(4,2),corners(4,1),'b+');
hold off
xlabel('Plot corner points found on briefcase');
end
%========End

%Now find the bounding box for this image
[verticesX, verticesY, centroidValue] = calcBoundingBox(openedImage);
bbWidth = verticesX(2) - verticesX(1);
bbHeight = verticesY(3) - verticesY(1);

%Find the vertices for the image projection
[rowVertices, colVertices] = find(openedImage == 1);

%Check for incorrect corner detections. If incorrect, then use the
%bounding box coordinates
if norm(corners(1,:) - centroid)>100
  corners(1,:) = [verticesY(1) verticesX(1)];
end
if norm(corners(2,:) - centroid)>100
   corners(2,:) = [verticesY(1) verticesX(2)];
end
if norm(corners(3,:) - centroid)>100
   corners(3,:) = [verticesY(3) verticesX(2)]; 
end
if norm(corners(4,:) - centroid)>100
   corners(4,:) = [verticesY(3) verticesX(1)]; 
end
%Plot points using the bounding box coordinates if
%necessary
if displayImages>0
% subplot(3,3,2);
imshow(openedImage);
hold on
plot(corners(1,2),corners(1,1),'b+');
plot(corners(2,2),corners(2,1),'b+');
plot(corners(3,2),corners(3,1),'b+');
plot(corners(4,2),corners(4,1),'b+');
hold off
xlabel('Plot corner points found on briefcase - After Bounding Box');
end

%Display the bounding box
if displayImages
%     subplot(3,3,3);
    imshow(openedImage);
    hold on 
    rectangle('Position', [verticesX(1), verticesY(1),...
              bbWidth, bbHeight], 'EdgeColor', 'red');
    hold off
end
% corners = 
%******************************************************************
%Perform an image transfer from the video frame
%to the original image
%*****************************************************************

% get input video frame and its size
inimage=videoFrame;
[IR,IC,D]=size(inimage);

%Get the portion of the video to place on the briefcase
UV=zeros(4,2);
XY=zeros(4,2);

%Need to find the target points manually
UV=corners;% target points
    
%The dimensions of the image. 
XY=[[1,1]',[1,IC]',[IR,IC]',[IR,1]']';    % source points

P=esthomog(UV,XY,4);    % estimate homography mapping UV to XY
v=zeros(3,1);
 %at each loop 
    for i=verticesY(1):verticesY(3)
        for j=verticesX(1):verticesX(2)
            point = xyzImage(i,j,:);
            point = permute(point,[3 1 2]);
            dist(i,j) = getEuclideanDistance(point,plane);
            
            v=P*[i,j,1]';        % project destination pixel into source
            y=round(v(1)/v(3));  % undo projective scaling and round to nearest integer
            x=round(v(2)/v(3));
            
            if dist(i,j)<0.1 && dist(i,j)~=-1 && (x >= 1) && (x <= IC) && (y >= 1) && (y <= IR)
                %overlay the color from the background image
                rgbImage(i,j,:) = inimage(y,x,:);
            end
        end
    end
    rgbimageFinal = rgbImage;
%******************************************************************


if displayImages>0
    figure(3)
   imshow(rgbImage);
   xlabel('The image with superimposed briefcase before removing noise');
end